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ABSTRACT 


Shape representation has been a central issue in computer 
graphics and computer-aided geometric design. Many physical 
phenomena involve curves and surfaces that are monotone (in 
some directions) or are convex. The corresponding representation 
problem is given some monotone or convex data, and a monotone 
or convex interpolant is found. Standard interpolants need not 
be monotone or convex even though they may match monotone or 
convex data. 

Several authors including Fritsch, Butland, McAllister, 
McLaughlin, Schumaker and Foley have investigated this pro- 
blem. Most of their methods involve the utilization of quad- 
ratic splines or Hermite polynomials. In our investigation, we 
have adopted a similar approach. These methods require derivative 
information at the given data points. The key to the problem is 
the selection of the derivative values to be assigned to the 
given data points. Schemes for choosing derivatives were exam- 
mined. Along the way, fitting given data points by a conic sect- 
ion has also been investigated as part of the effort to study 
shape-preserving quadratic splines. 



I. INTRODUCTION 


This project commenced in May 1987 and ended in October 1988, 
which includes a six-month, no-cost extension. Dr. Maria Lam 
was the Principal Investigator. She was assisted by a graduate 
student in the Applied Mathematics program at Hampton University. 

Hampton University has benefited by this project. Through 
this grant, the Principal Investigator was able to obtain release 
time and physical resources to conduct research into the area 
that interests her. This project has been beneficial to the 
graduate program in Applied Mathematics by way of the research 
assistantship made available to the graduate student. The work 
done under this project will constitute part of the graduate 
student's dissertation. The Computer Science Department has also 
benefited. This project helped to fulfil part of the faculty 
research criterion set forth by the ACM Accreditation Board when 
the Department applied for accreditation in the spring of 1988. 

In short, this project has been beneficial to the University in 
several ways. 

Two contributed papers were presented on research done under 
this project. The details of these are given in Sections II and 


III. 


II. FITTING DATA BY CONIC SECTIONS 


Data Identification : 

The initial attempt has been to identify 2-D data assuming 
there is some underlying relation among the data points. Owing 
to the diversity of this problem, the study has focused on conic 
sections. A conic section can be represented either algebraically 
or geometrically /l,2/. Some of these representations require 
additional information than just the data points. In the absence 
of this information, in theory, a conic section can be represent- 
ed algebraically by a quadratic relation 

ax 2 +2bxy+cy 2 +2dx+2ey+f=0 ( 1 ) 

or PQP t =0 where 

P=[x,y,l] 

a b d 
Q= bee 
. d e f _ 

Define angle s as 

2s=tan _1 (2b/ (a-c) ) 
or 

s-TC/4 if a=c. 

Now effect a rotation of the coordinate system given by 
x=Xcos(s) - Ysin(s) 
y=Xsin(s) + Ycos(s). 

Thus, equation (1) becomes 
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AX 2 +CY 2 +2DX+2EY+F=0 


( 2 ) 


By examining the coefficiets of equation ( 2 ), the type of conic 
section (including the degenerate cases) can be determined /3/. 
Hence a conic section is completely determined by five data 
points since only five of the six coefficients in equation (1) 
are independent . 

One can indirectly solve for these coefficients by the foll- 
owing method /4/. Assume the five data points are (x^y^) , 
( x i+l»yi+l) • ( x i+ 2 /Yi+ 2 )' ( x i+ 3 'Yi+ 3 ) and ( x i+ 4 'Yi+ 4 )‘ Form the 

6x6 matrix M given by 


X 1 

x iYi 

Yi 

x i 

Yi 

1 ' 

x l+l 

x i+lYi+l 

Yi+1 

x i+i 

Yi+1 

1 

x l+2 

x i+2Yi+2 

Yl+2 

x i+2 

Yi+2 

1 

x l+3 

x i+3Yi+3 

yf+3 

x i+3 

Yi+3 

1 

x l+4 

x i+4Yi+4 

Yi+4 

x i+4 

Yi+4 

1 

X 2 

X Y 

y 2 

X 

y 

1 _ 


The equation 

det(M)=0 (3) 

yields a quadratic equation in x, y. When x=xj c and y=yj c 
(i<=k<=i+4) , two of the rows of the matrix are identical, and 
equation (3) is satisfied. This implies that equation (3) 
represents the unique conic passing through these five points. 

However, this algebraic method is found to be numerically 
unstable. The determination of the type of conic section depends 
largely on whether certain coefficients are positive, negative or 
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zero. It is extremely sensitive to small perturbations in coeff- 
icients, which may alter the type of the conic section. When it 
is implemented in a computer program, the situation is even 
worse. Such method routinely involves considerable amount of 
algebraic operations. Their floating-point implementation on a 
computer almost invariably produces perturbations in the coeff- 
icients. The selection of reasonable tolerance is difficult. 
Computer algebras such as MACSYMA and SMP are ideal tools for 
this type of investigations because actual evaluation of 
expressions is delayed until the final step. This leads to 
improved accuracy, but, however it does not eliminate the 
problem. 

Several sets of conic data are generated. For each data set, 
five consecutive data points are used for fitting. A rotation is 
performed to eliminate the mixed term. Then the ceofficients of 
the result are observed. This process is repeated until all 
points of a data set are exhausted. 

Results : 

The results of this study are mixed. As one would expect, it 
is difficult to identify these figures definitely. 

To be identified as data that come from a circle, the gener 
ated coefficients A, C must be equal. This rarely happens. The 
best one can expect is to identify the data as elliptic. Our 
study shows that one set of circular data is identified as ell ip 
tic. A second set is identified as elliptic and hyperbolic with 
the elliptic outnumbering the hyherbolic. A third data set prov- 
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ides similar result. Elliptical data exhibit similar behavior. 

Parabolas are identified by A=0 or C=0 but not both. For one 
set of parabolic data with the values of x several units apart, 
we are able to identify this trend every time. A second set of 
data, with closely spaced x values, does not fare well at all. 

It is correctly identified as parabola only on two occasions. For 
most part, it has been identified incorrectly as hyperbolic 37 
times and elliptic 5 times. The third set is the most curious 
one. The first 18 fits indicate the data is hyperbolic, this is 
followed by 22 fits that indicate poarabolic data. 

If the data is hyperbolic, the result of algebraic method 
leads one to conclude that the data is either hyperbolic or 
elliptic; however the former conclusion (i.e. hyperbolic) occurs 
more frequently. 

In spite of this lack of positive identification, one pattern 
emerges. Parabolic data are identified as either parabolic or 
hyperbolic. It is the only type of data that will draw parabolic 
conclusion. Circles and ellipses are classified as ellipses or 
hyperbolas with the elliptical conclusion appearing more often 
than the other. This allows one to draw rough estimates on the 
type of conic section that the data could have originated from. 

In summary, when given a set of conic data and algebraic 
method is used in attempt to identify it, one would assume the 
data is hyperbolic initially. If on some occasions the data is 
identified as parabolic and as hyperbolic on other occasions, 
one should accept the data as parabolic. If the data is ident- 
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ified as ellipse and hyperbola with ellipse more frequent than 
hyperbola, one should conclude it to be elliptic. 

Data Fitting : 

Our next attempt is to fit arbitrarily given data sets by 
conic section. Several data sets including some standard data 
sets to test shape-preserving interpolants are used. The proced- 
ure described in the previous subsection minus the rotation is 
repeated on these data sets. At the end of each fitting, the 
derivative at the point P(i+2) is calculated. These derivative 
values will be examined later. 

When conic data is used, the result is rather consistent if 
consecutive data points are no more than a few units apart. The 
equations of the fits are approximately the same. The changes 
in coefficients are consistent and gradual. However, we 
are unable to recover the equation used to generate the data. 
These fits do not resemble the original equation whatsoever. For 
other data sets, if the five points used to do the fitting do not 
conform to conic shape, the fit will be poor if not downright 
impossible. This is usually indicated by some coefficients having 
values close to zero. The derivatives obtained from these conic 
fits generally preserve the monotone property of the data set. 
They also preserve the convexity of nice, well-behaved data sets, 
but fail under severe testing when this method is applied to some 
standard data sets used to test shape-preserving property of 
interpolants. 


6 


Data Used : 

The following are some of the data sets used in this invest 
igation. 

Set 1 


(y-3) 2 = -2 Ox 


X -79 

y 42.75 

-76 

41.99 

-73 -72 -70 -66 -63 

41.21 41 40.42 39.33 38.5 

-62 

38.21 

Set 2 




x 2 +y 2 = 100 2 




x 84 

100 
y 54 

0 

86.6 

50 

89.1 93.4 95.1 97.8 

45.4 35.8 30.9 20.8 

99.5 

9.99 


Set 3 

x 2 +y 2 = 100 2 

X 84 85.4 

86.6 

88 

89.1 

90 

92 

93.4 94 

95.1 

96 

97 

97.8 

98.5 

99 99.5 

100 





y 54 52 

50 

47.5 

45.4 

43.59 

39.19 

35.8 34.12 

30.9 

28 

24 .31 

20.8 

17.25 

14.11 9.99 

0 






Set 4 

y = x 3 -x+l 

X 

0 

1 

2 

2.25 

2.8 

3 

3.1 


3.6 

4 






y 

i 

1 

7 

10.14 

20.15 

25 

27.69 


44.06 

61 
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set 5 

McAllister and Roulier data 


x 0 

2 

4 

6 

8 

10 

12 

y o 

2 

44 

88 

132.1 

1132.1 

2132.2 

Set 6 

Roulier data 








XO 1 2 3 4 5 6 

7 

y 0 1.8 2.7 2.9 3 3 3.65 

4.7 


Set 7 

Pruess data 


X 

22 

22.5 

22.6 

22.7 

22.8 

22.9 

23 


23.1 

23.2 

23.3 

23.4 

23.5 

24 


y 

523 

543 

550 

557 

565 

575 

590 


620 

860 

915 

944 

958 

986 



Set 8 

Akima data 

x 0 

2 

3 

5 

6 

8 

9 

11 

12 

14 

15 




y 10 

10 

10 

10 

10 

10 

10.5 

15 

50 

60 

85 





Set 

RNP 

9 

14 data 

(Fritsch and Carlson) 



X 

7.99 

8.09 

8.19 

8.7 


9.2 

20 

10 

12 

15 

y 

0 

0.0000276429 

0.0437498 

0.169183 


0.469428 

0.999994 

0.94374 

0.998636 

0.999919 
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III. SHAPE-PRESERVING INTERPOLATION 


The problem is: Given points x 1 <x 2 <...<x n and values {yj.}i=i, 
find f such that 

f ( x i) = Yi i=l ,2, . . . ,n 

and f preserves the monotonicity and/or convexity of the data, f 
is usually found by piecing fj/s together where f^ is defined 
over [Xi,Xi +1 ] and 

fi(*i) = Yi, fi(*i + i) = Yi+i- 

Let h x x , d (y y j_) /h j_ , i=l,2, . . . ,n— 1. 

There are two general approaches to this problem. One employs 
geometric constructions which could be quite involved /5,6/. The 
other approach uses quadratic or cubic splines, together with 
carefully chosen derivative values or tension parameters /7,8,9/ 
to control the shape of the curve. Quadratic spline is chosen for 
this investigation because quadratic polynomials have one less 
extremum than cubic polynomials, hence better preserve the 
geometric properties of the given data. The derivatives of 
quadratics are linear, this makes them easier to study and 
manipulate. Hence we focus on a simple two-point Hermite inter- 
polation problem involving quadratic polynomials: to find a 
quadratic polynomial h, if possible, such that 

h(*i)=yi, h'(x i )=s i , i-1,2 (4) 

for some preassigned values 
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Conditions on the Derivatives : 


Schumaker /6/ proved that there exists a quadratic poly- 
nomial solving this problem if and only if 

s l +s 2 _ Y2-Y1 

2 3C 2 ~x 1 


(5) 


and explicitly defined this polynomial. Moreover, this poly- 
nomial is monotone on I=[Xi,X 2 ] if SiS 2 >= 0 , it is convex on 
I if s 1 <s 2 / and concave if s 1 >s 2 . 

When condition (5) fails, Schumaker showed that corres- 
ponding to every u in (X;l,x 2 ), it is possible to choose a 
functional value y* and derivative s* at u, and define a 
quadratic spline over I that satisfies condition (4) . If 
s l s 2 >= 0, then it is monotone on I if and only if sis*>=0. 

If s 1 <s 2 , it is convex on I if and only if s 1 <=s*<=s 2 . 
Similarly, if s 1 >s 2 , then it is concave if and only if 
s 1 >=s*>=S2 . 


Selection of Derivatives : 

The conditions stated in the previous subsection are 
deceptively simple. However, when both monotonicity and 
convexity are taken into consideration, the problam takes 
on a new dimension because one cannot satisfy one condition 
at the expense of the other. One generally approaches this 
problem by first assigning some carefully chosen derivative 
value to each data point, then modify some of these values 
if necessary to satisfy the monotone and/or convexity 
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conditions, then the knots are inserted and their functional 
values and derivatives are determined, and finally the coeff- 
icients of the quadratic polynomials constituting the spline 
are calculated. 

Schumaker /6/ used a weighted average of d^_^ and d^ 
as the value of s-^. McAllister and Roulier /5/ used a very 
involved geometric construction to obtain their shape- 
preserving quadratic splines. But their method can be 
identified as Schumaker' s formulation with start-up deri- 
vative s^ defined as the harmonic mean of di_ 1# d^ and a 
modification scheme for s^. DeVore and Yan /8/ improved 
upon their modification. 

This investigation has focused on the selection of 
the values of s^. The following method looks promising. 

The derivative s^ is defined as: 


di-idi 

Si = < tdi_ 1 +(l-t)di 

l 0 


where 


if di_idi>0 
otherwise 


t 


1 


- ( 1 + 
3 


h i_l+hi 


). 


This formula is symmetric in d^_ lf d^ because one can simplify 
(1-t) to obtain 

1 hi 

1-t =1 ( 1 + ) 

3 h i-i +h i 


li 


1 



3 

1 

3 h-£_2. + Ii^ 

1 hi.! 

- - ( 1 + ) • 

3 h i-l+hi 

This makes the formula independent of the order of the data 
points, i.e. the formula gives the same value when the data 
is processed from the left or from the right. This Si also 
has the properties: 

1. minfdi.i^i) <=si<=max(di_i,di) 

thus the slope of the curve lies between the slopes of 
the two adjacent data segments, 

2. | Si| <=min(3 | di.i | , 3 |di | ) and 

3. Si reduces to the harmonic mean of di.^, di when xi_i, 

x i' x i+l are equally spaced (i.e. when hi_i=hi) . 

All these properties make it very desirable. 

The investigation was still in progress when this project 
came to an end. It will be carried on under the new project 
as grant NAG-1-948. 
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IV. CONTRIBUTED PAPERS 


The following papers were presented on the research 
conducted under this project. 

1. Shape Preserving Interpolants 
Maria H. Lam 

Annual Meeting of Virginia Academy of Science, Norfolk, 
Virginia, May 1987 

Virg. J. Sc., Vol . 38, No. 2, 71 (1987) 

2 . Fitting Data by a Conic Section 
Maria H. Lam 

Annual Meeting of Virginia Academy of Science, 
Charlottesville, Virginia, May 1988 
Virg. J. Sc., Vol. 39, No. 2, 103 (1988) 

Abstracts of these papers are shown on following pages. 
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SHAPE PRESERVING INTERPOLANTS 
Maria H. Lam 

Department of Computer Science 
Hampton University 
Hampton, VA 23668 

Recently there has been immense interest in shape-preserving 
interpolations. Most of these algorithms are quite involved. H. 
McLaughlin has proposed a simple, local method to generate inter 
polants which preserve the monotonicity and convexity of the 
given data. F. Fritsch and J. Butland have proposed a method to 
select values of derivatives at data points such that local mono 
tone piecewise cubic interpolants will be obtained. These two 
methods are implmented by using computer algebra SMP. They are 
applied to several sets of data. Curves produced by these two 
methods are compared. Some of their properties are discussed. 
(Supported by NASA under grant NAG-1-415 and NAG-1-760) 
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FITTING DATA BY A CONIC SECTION 


Maria H. Lam 

Department of Computer Science 
Hampton University 
Hampton, VA 23668 


ABSTRACT 


A conic section can be represented either algebraically or 
geometrically. In theory, it can be represented algebraically 
by a quadratic relation of the form 

ax 2 +2bxy+cy 2 +2dx+2ey+f=0 . 

It is completely determined by five data points as only five of 
the six coefficients are independent. These coefficients can be 
determined by solving a system of five linear eguations in six 
unknowns indirectly. This method is found to be numerically 
unstable. Its floating point implementation on a computer almost 
invariably produces perturbations in the coefficients which may 
alter the character of the conic section. Several sets of conic 
data are generated. For each data set, five consecutive data 
points are used for fitting. A rotation is performed to eliminate 
the mixed term. The coefficients of the result are observed and 
compared. This process is repeated until all points of a data 
set is exhausted. Result of this research will be presented. 
(Supported by NASA under grant NAG-1-760) 
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